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LONG-TERM  GOALS 

•  To  develop  a  method  of  modeling  sound  propagation  in  an  environment  with  multi-scale 
inhomogeneities,  which  preserves  the  efficiency  and  intuitive  qualities  of  the  ray  theory  but  is 
free  from  spurious  environmental  sensitivity  and  strong  perturbations  associated  with  ray 
trajectories. 

•  To  investigate  and  quantify  effects  on  underwater  acoustic  wavefronts  of  internal  gravity 
waves,  sea  swell,  “spice,”  and  other  small-scale  processes  in  the  water  column. 

OBJECTIVES 

1 .  To  assess  significance  of  time  dependence  of  the  sound  speed  and  flow  velocity  perturbations 
on  predictability  of  acoustic  wavefronts  and  timefronts. 

2.  To  quantify  horizontal  refraction  of  sound  by  random  meso-scale  inhomogeneities  at  0(l)Mm 
propagation  ranges. 

3.  To  find  the  variance  and  bias  of  random  ray  travel  times  in  the  regime,  where  the  ray 
displacement  may  be  comparable  to  the  vertical  extent  of  the  underwater  waveguide  but  the 
clustering  has  not  developed  yet. 

4.  To  determine,  using  a  perturbation  theory  and  numerical  simulations,  typical  propagation 
ranges  where  clustering  of  chaotic  rays  replaces  the  anisotropy  of  ray  scattering  as  the  main 
physical  mechanism  responsible  for  acoustic  wavefront  stability. 

5.  To  develop  an  efficient  technique  for  modeling  acoustic  wavefronts  and  their  perturbations  in 
range-dependent  and  horizontally  inhomogeneous  oceans. 

6.  To  model  perturbations  of  acoustic  wavefronts  and  timefronts  by  internal  gravity  waves, 
internal  tides,  sea  swell,  and  “spice”  in  the  ocean. 

7.  To  investigate  implications  of  wavefront  stability  on  the  downward  extension  of  acoustic 
timefronts  and  deepening  of  lower  turning  points  of  steep  rays  due  to  small-  and  meso-scale 
physical  processes  in  the  upper  ocean. 
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APPROACH 


Our  primary  theoretical  approach  is  an  extension  of  the  method  employed  in  (Godin,  2003,  2007, 
2009),  where  a  novel  perturbation  technique  has  been  developed  to  solve  the  eikonal  equation  and 
calculate  wavefront  and  ray  trajectory  displacements,  which  are  required  to  be  small  over  a  correlation 
length  of  the  environmental  inhomogeneities  but  not  necessarily  over  the  entire  acoustic  propagation 
path.  Methods  based  on  the  geometrical  acoustics  have  been  supplemented  by  full-wave  approaches 
(Brekhovskikh  and  Godin,  1999),  including  a  normal-mode  theory  for  range-dependent  and 
horizontally  inhomogeneous  waveguides. 

The  analytical  methods  have  been  complemented  by  numerical  wavefront-tracing  techniques.  A  major 
problem  with  a  direct  modeling  of  acoustic  wavefronts  in  the  ocean  through  numerical  solution  of  the 
eikonal  equation  lies  in  the  eikonal  (and  acoustic  travel  time)  being  a  multi-valued  function  of 
position.  A  number  of  computational  approaches  to  solve  the  eikonal  equation  without  ray  tracing  had 
been  developed  in  mathematical  and  seismological  communities  (Vidale,  1990;  Sava  and  Fomel, 

2001;  Sethian,  2001,  Benamou,  2003).  However,  most  of  these  methods  are  only  capable  of  tracing  the 
wavefront  of  the  earliest,  first  arrival  and  thus  are  not  suitable  for  underwater  acoustic  applications. 

We  have  adapted  to  ocean  acoustics  problems  the  Lagrangian  wavefront  construction  techniques 
(Vinje  et  ah,  1993,  1999;  Lambare  et  ah,  1996;  Sava  and  Fomel,  2001;  Chambers  and  Kendall,  2008; 
Hauser  et  ah,  2008),  which  were  originally  developed  in  the  context  of  exploration  seismology  and  are 
capable  of  computing  all  arrivals.  Extension  of  the  existing  wavefront  construction  techniques  to  long- 
range  underwater  sound  propagation  is  a  non-trivial  and,  in  fact,  rather  challenging  task  because  the 
number  of  ray  arrivals  and  topological  complexity  of  wavefronts  in  the  ocean  far  exceed  those  in  the 
geophysical  applications  considered  to  date. 

Theoretical  results  and  new  modeling  capabilities  have  been  verified  against  numerical  simulations 
perfonned  with  well-established  ray  and  PE  propagation  codes. 

The  key  individuals  that  have  been  involved  in  this  work  are  Oleg  A.  Godin  (CIRES/Univ.  of  Colorado 
and  NOAA/ESRL),  Nikolay  A.  Zabotin  (CIRES/Univ.  of  Colorado),  and  L.  Y.  Zabotina  (ECEE/Univ. 
of  Colorado).  Dr.  Zabotin  and  Ms.  Zabotina  focused  on  developing  and  testing  an  efficient  computer 
code  for  modeling  acoustic  wavefronts  in  inhomogeneous,  moving  fluids.  Dr.  Godin  took  the  lead  in 
theoretical  description  of  effects  of  the  ocean  currents,  localized  inhomogeneities,  and  sound-speed 
time -dependence  on  the  acoustic  wavefronts. 

WORK  COMPLETED 

Most  of  the  work  on  the  project  was  completed  by  the  end  of  FY’  13.  In  FY’  14,  our  3-D  Huygens 
Wavefront  Tracing  (3-D  HWT)  algorithm  was  applied  to  study  effects  of  various  oceanographic 
processes  on  underwater  acoustic  fields.  The  3-D  HWT  algorithm  was  also  applied  to  investigate  long- 
range  propagation  of  infrasound  in  the  atmosphere  (Zabotin  et  ah,  2014)  and  of  infragravity  waves  in 
the  ocean  (Godin  et  ah,  2014).  Infragravity  waves  are  surface  gravity  waves  with  periods  from  ~0.5  to 
-50  minutes.  Infragravity  waves  are  refracted  by  ocean  currents  and,  especially,  by  seafloor 
bathymetry  and  can  be  treated  in  horizontally  inhomogeneous  ocean  using  the  vertical  modes  - 
horizontal  rays  theory  [see,  e.g.,  (Brekhovskikh  and  Godin,  1999)].  The  3-D  HWT  algorithm  was  used 
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to  simulate  horizontal  structure  of  infragravity  wave  fields  and  interpret  observations  of  seafloor 
pressure  fluctuations  (Godin  et  ah,  2014). 

RESULTS 

A  software  package  was  developed  for  numerical  modeling  of  3-D  sound  propagation  in  the  ocean, 
where  sound  speed  and  current  velocity  are  arbitrary  smooth  functions  of  three  spatial  coordinates.  The 
package  is  freely  available  through  the  Ocean  Acoustics  Library  at 

http://oalib.hlsreseareh.eom/Rays/HWT_3D_mm_v.l.0/  .  The  software  package  is  a  numerical 
implementation  of  an  alternative  formulation  of  the  geometrical  acoustics,  where  wavefronts  rather 
than  rays  are  the  geometric  “backbone”  of  the  wave  field.  Wavefront  tracing  can  be  viewed  as  a  direct 
implementation  of  the  Huygens’  principle.  Our  current  wavefront  tracing  algorithm  (Zabotin  et  ah, 
2014)  is  an  extension  to  three-dimensionally  inhomogeneous,  moving  media  of  an  earlier  2-D 
algorithm  for  sound  in  moving  fluids  (Zabotin  et  ah,  2012)  and  the  original  Huygens  Wavefront 
Tracing  (HWT)  algorithm  (Sava  and  Fomel,  2001),  developed  as  a  part  of  the  open-source  Madagascar 
project  for  seismic  modeling  and  imaging.  Aside  from  generalization  to  moving  media,  the  extended 
HWT  algorithm  improves  the  accuracy  and  stability  of  the  original  algorithm  and  includes  new 
functionality  to  account  for  sound  reflections  at  sloping  boundaries  and  to  efficiently  model  acoustic 
timefronts  in  addition  to  the  wavefronts.  For  wavefront  tracing  in  inhomogeneous  media,  HWT  is 
much  more  computationally  efficient  and  robust  than  traditional  ray  codes  (Sava  and  Fomel,  2001). 
Unlike  many  other  eikonal  solvers,  the  HWT  method  produces  the  output  in  ray  coordinates  and  has 
the  important  ability  to  track  multiple  arrivals.  With  the  HWT  method,  each  wavefront  is  generated 
from  the  preceding  one  by  finite  differences  in  the  ray-coordinates  domain. 

The  3-D  HWT  algorithm  has  been  applied  to  study  effects  of  various  oceanographic  processes  on 
underwater  sound  propagation  and  also  has  been  demonstrated  to  be  an  efficient  and  robust  technique 
for  modeling  infrasound  propagation  in  the  atmosphere  (Zabotin  et  ah,  2014). 

To  illustrate  application  of  the  3-D  wavefront  tracing  to  modeling  sound  propagation  in  shallow  water,  we 
consider  environmental  conditions  similar  to  that  of  the  1995  Shallow-Water  Acoustic  Random  Media 
(SWARM)  experiment  (Apel  et  ah,  1997).  Rather  strong  sound  speed  variations  in  time  and  in  the 
horizontal  plane  were  encountered  during  the  SWARM  experiment.  The  sound-speed  variability  was 
primarily  due  to  packets  of  nonlinear  internal  gravity  waves  (IWs),  which  were  generated  by  tides  on  a 
shelfbreak.  We  follow  (Apel,  2003;  Godin  et  ah,  2006)  and  use  the  “dnoidal”  model  of  the  vertical 
isopycnal  displacement  in  the  nonlinear  IW  packets.  In  our  simulations  we  placed  the  source  of  sound  at 
the  point  with  coordinates  x  =  0,  y  =  5000  m,  z  =  20  m.  The  source  depth  corresponds  to  the  peak  of  the 
first  IW  mode  shape  function.  Distortions  of  a  particular  wavefront,  which  corresponds  to  the  acoustic 
travel  time  of  6  s,  are  shown  in  Figure  1  for  different  positions,  i.e.,  stages  of  evolution,  of  the  IW  packet. 
The  IW  packet  positions  are  characterized  by  distance  D  =  Vt  the  packet  covered  since  its  generation;  D  = 
0  corresponds  to  sound  propagation  in  the  absence  of  IWs,  while  D  =  30  km  corresponds  to  the  time,  when 
the  peak  of  the  sound  speed  distortion  reaches  the  y  =  5000  m  plane.  Note  that  the  acoustic  wavefront 
maintains  its  general  shape  during  the  IW  packet  propagation  despite  strong  perturbations  of  individual 
rays  (not  shown)  that  form  the  wavefront. 

Relative  variations  of  sound  speed  in  the  atmosphere  are  much  larger  than  in  the  ocean.  Unlike  the 
ocean,  where  the  Mach  number  of  currents,  understood  as  the  ratio  of  the  flow  velocity  to  the  sound 
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speed,  rarely  exceeds  10  ,  winds  can  have  Mach  numbers  as  large  as  0.5  and  are  a  major  factor 
shaping  acoustic  fields.  Combined  with  strong  variations  of  the  sound  speed  and  wind  velocity  with 
height,  this  leads  to  a  rather  complicated  arrival  structure  of  acoustic  fields  due  to  distant  compact 
sources.  Understanding  and  utilizing  acoustic  arrivals  in  the  atmosphere  requires  an  accurate  modeling 
of  sound  propagation  in  moving  fluids.  We  illustrate  atmospheric  applications  of  the  HWT  algorithm 
by  modeling  propagation  of  infrasound  generated  by  Eyjafjallajokull  volcano  in  southern  Iceland. 
Eruptions  of  this  volcano  were  recorded  by  infrasonic  sensors  at  distances  up  to  3700  km  in  2010 
(Drob  et  ah,  2010). 

Modern  atmospheric  models,  such  as  the  Whole  Atmosphere  Community  Climate  Model  (WAM)  (Liu 
et  ah,  2010;  Akmaev,  2011),  provide  information  about  three-dimensional  distributions  of  all  major 
atmospheric  parameters  necessary  for  sound-speed  calculations:  the  temperature  T,  zonal  (w)  and 
meridional  (v)  components  of  the  wind  velocity  as  well  as  the  content  of  water  vapor  and  other 
chemical  constituencies.  To  specify  the  propagation  medium,  we  used  a  version  of  WAM  output  that 
provides  ~2  deg  resolution  in  geographic  coordinates  and  150  unevenly  distributed  grid  points  in  the 
altitude,  from  0  to  500  km  above  sea  level  (Figure  2).  Figures  2b,  c  illustrates  the  complicated  spatial 
structure  of  the  horizontal  and  vertical  variations  of  the  physical  parameters  of  the  atmosphere  that 
govern  infrasound  propagation  in  the  Eyjafjallajokull  volcano  vicinity. 


8850  8900  8950  9000  9050  9100 
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Figure  1.  Evolution  of  an  acoustic  wavefront  during  propagation  of  a  packet  of  nonlinear  internal 
waves.  Cross  section  by  the  vertical  plane  y  =  5  km  of  the  acoustic  wavefront  at  travel  time  t  =  6  s  is 
shown  for  five  different  stages  of  the  internal  wave  packet  development.  A  point  sound  source  is 
located  at  (0,  5000  m,  20  m).  The  internal  wave  train  originates  at  a  location  with  horizontal 
coordinates  x  =  5  km,  y  =  35  km.  The  packet  occupies  a  circle  with  radius  D  =  0,  29  km,  30  km,  or 
32  km.  When  D  =  0,  sound  propagation  is  unaffected  by  internal  waves.  When  D  =  30  km,  the 
maximum  of  the  sound  speed  perturbations  reaches  they  =  5  km  plane. 

[The  acoustic  wavefront  consists  of  a  number  of  connected  branches,  most  of  which 
span  the  entire  waveguide  from  the  seafloor  to  the  ocean  surface.  The  wavefront 
maintains  its  general  shape  during  IW packet  propagation.  It  acquires  additional 
branches  and  extends  to  larger  x  values  for  smaller  D  values.] 
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Figure  2.  Atmospheric  model  used  in  simulations  of  infrasound  propagation,  (a)  Global  distribution 
of  the  absolute  air  temperature  (in  Kelvin)  at  altitude  10  km  predicted  by  the  Whole  Atmosphere 
Community  Climate  Model  (WAM).  Eyjafjallajokull,  Iceland,  location  is  marked  with  an  asterisk, 
(b)  Vertical  profdes  of  the  temperature  T,  sound  speed  c,  zonal  u  and  meridional  v  components  of 
wind  velocity  at  the  volcano  location,  (c)  Vertical  profdes  of  the  Mach  number  of 

wind  at  the  volcano  location. 

[Globally,  absolute  air  temperature  varies  between  195  K  and  240  K,  with  the  coldest  regions 
between  75  and  90  degrees  southern  latitude  and  the  hottest  regions  at  15-35  degrees  of  southern 
and  northern  latitudes  in  the  eastern  hemisphere.  Above  the  volcano,  the  temperature  and  sound 
speed  vary  rapidly  below  100  km  and  steadily  increase  with  height  above  100  km.  Magnitudes  of  the 
zonal  and  meridional  components  of  the  wind  velocity  reach  130  m/s  and  290  m/s,  respectively.  The 
maximum  value  of  the  Mach  number  of  the  wind  is  about  0.3  below  250  km  and  0.5  between  250 

and  500  km.] 


In  acoustic  simulations,  we  disregarded  topographic  effects  and  curvature  of  the  Earth  and  modeled 
ocean  and  ground  surface  as  a  horizontal  plane.  With  the  origin  of  coordinates  chosen  at  sea  level  at 
the  location  of  the  Eyjafjallajokull,  horizontal  coordinate  axes  Ox  and  Oy  were  directed  towards  east 
and  north,  respectively.  Although  not  resolved  on  the  scale  of  the  figures,  almost  all  ray  paths  in 
Figure  3a  and  wavefronts  in  Figures  3b  and  3  c  are  duplicated  because  of  the  reflection  from  the 
ground  in  the  vicinity  of  the  elevated  sound  source.  Most  of  the  other  peculiarities  of  the  ray  and 
wavefront  geometry  are  due  to  irregular  three-dimensional  wind  distribution  predicted  by  WAM.  Note 
the  very  strong  anisotropy  of  infrasound  propagation  in  the  horizontal  plane,  which  can  be  appreciated, 
in  particular,  from  the  shape  of  a  wavefront  on  the  ground  (Figure  3  c). 

The  critical  role  of  winds  in  infrasound  propagation  demands  a  careful  mathematical  description  of  the 
effects  of  the  fluid  flow  on  acoustic  fields.  Effective  sound  speed  approximation,  where  a  moving 
medium  is  replaced  by  a  motionless  fluid  with  an  effective  sound  speed,  is  widely  used  in  ocean  and 
atmospheric  acoustics.  Our  HWT  algorithm  does  not  rely  on  the  effective  sound  speed  approximation 
and  allows  one  to  quantify  its  accuracy. 
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Figure  3.  Application  of  the  3-D  HWT  algorithm  to  numerical  modeling  of  infrasound  propagation 
in  the  atmosphere,  (a)  Projection  on  the  vertical  plane  y  =  0  of  rays  leaving  a  point  source  at  the 
azimuthal  angle  i//  =  0  and  with  different  grazing  angles.  The  point  source  is  placed  at  the  summit  of 
the  Eyjafjallajokull  volcano  at  elevation  Zo  =  1666  m  above  sea  level,  (b)  Cross  section  of  the 
acoustic  wavefront  at  t  =  800  s  after  an  eruption  by  the  vertical  plane  y  =  0  passing  through  the 
source.  In  addition  to  the  acoustic  wavefront  in  the  atmosphere  considered  as  a  moving  fluid,  the 
result  of  calculation  of  the  same  wavefront  in  an  effective ”  motionless  fluid  is  shown,  (c) 
Intersections  of  the  wavefronts  shown  in  (b)  with  sea  and  ground  surfaces  modeled  as  the  horizontal 
plane  z  =  0.  For  comparison  of  the  accurate  and  approximate  wavefront  calculations,  the  insert 

shows  a  blow-up  of  a  portion  of  the  wavefronts. 

[A  projection  of  a  set  of  rays  launched  from  a  point  source  shows  that  rays  with  smaller  grazing 
angles  have  turning  points  and  can  return  to  the  ground  to  the  east  and  to  the  west  of  the  source. 
Turning  points  occur  up  to  heights  as  large  as  220  km.  Steep  rays  do  not  have  turning  points  and  do 
not  return  to  the  ground.  The  wavefronts  calculated  for  actual  and  “effective”  motionless  fluid 
become  drastically  different  at  high  elevations.  At  ground  level,  the  accurate  and  approximate 
calculations  give  similar  results  but  with  ~7  km  shift  between  the  actual  and  approximate 

wavefronts.] 


In  the  Eyjafjallajokull  eruption  scenario  we  consider,  the  most  striking  differences  between  wavefronts 
in  moving  and  “effective”  media  are  observed  at  large  elevations  (Figure  3b).  Unlike  the  actual, 
moving  medium,  a  large  shadow  zone  arises  at  near-zenith  directions  in  the  “effective”  medium.  The 
shadow  zone  leads  to  termination  of  the  cross  section  of  the  approximate  wavefront  at  large  heights 
(Figure  3b).  Farge  distortions  of  the  wavefront  position  occur  at  intermediate  heights.  As  expected 
from  theoretical  analysis  of  the  effective  sound  speed  approximation  accuracy  (Godin,  2002),  the 
approximation  performs  much  better  for  the  shallower  rays,  which  are  responsible  for  infrasound 
arrivals  on  the  ground  at  t  =  800  s  (Figure  3c).  For  the  portions  of  the  wavefront  shown  in  the  insert  in 
Figure  3c,  the  accurate  and  approximate  wavefronts  are  separated,  on  average,  by  about  1  km,  which 
corresponds  to  an  early  arrival  bias  (i.  e.,  an  underestimation  of  the  acoustic  travel  time)  by  ~  3  s  in  the 
effective  sound  speed  approximation. 

IMPACT/APPLICATIONS 

Wavefront  tracing  has  been  established  as  a  new  approach  to  modeling  sound  fields  in  deep  and 
shallow  water,  which  readily  accounts  for  3-D  propagation  effects  and  current-induced  acoustic 
anisotropy.  In  addition  to  modeling  acoustic  manifestations  of  internal  gravity  waves,  internal  tides, 
and  swell  in  the  ocean,  wavefront  tracing  provides  an  efficient  technique  for  simulating  long-range 
propagation  of  infrasound  in  the  atmosphere  and  infragravity  waves  in  the  ocean. 

RELATED  PROJECTS 

None. 
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